#ifndef LINEARSYSTEMSOLVER_H
#define LINEARSYSTEMSOLVER_H

#define IMKL

#include "ParticleSystem.h"
#include "ConstraintSystem.h"
#include <math.h>
#include <stdio.h>

#ifdef IMKL
#include <mkl_lapacke.h>  
#include <mkl.h>  

#pragma comment(lib,"mkl_lapack95_lp64.lib")
#pragma comment(lib,"mkl_rt.lib")

typedef struct  
{
	double *A;
	double *superb;
	double *s;
	double *u;
	double *v;
	double **cache;
	double **pseInv;
	lapack_int m;
	lapack_int n;
	lapack_int ldu;
	lapack_int ldvt;
	lapack_int lda; 
}SVD_param;
#endif

enum SOLVER{
	CONJUGATEGRADIENT,
	GAUSSJORDAN,
	SINGLEVALUEDECOMPOSITION
};


class LinearSystemSolver
{
public:
	LinearSystemSolver(void);
	LinearSystemSolver(ParticleSystem *ps, ConstraintSystem *cs);
	~LinearSystemSolver(void);
	void Clear();
	void SystemInitialize(ParticleSystem *ps, ConstraintSystem *cs);
	void SystemConstruction();
	void MassBlockConstruction();
	void CGradBlockConstruction();
	void CGradTransposeBlockConstruction();
	void BVectorConstruction();
	void ExtForceConstruction();
	void HarmonicOscilatorConstruction();
	void Solve();
	void SystemIterativeSolve();
	void JacobiIter();
	void ConjugateGradIter();
	void GaussJordan();
	void SVDSolver();
	void PseudoInverse();
	void SVD();
	void GetAccleration(double *accleration);
	inline int I(int i, int j);

public:
	ParticleSystem *ps;
	ConstraintSystem *cs;
	double **A;
	double **CGrad;
	double *b;
	double *residual;
	double *cache;
	double *residual_cache;
	double *p;
	double *HO;
	double *output;
	double **GS;
	int scale;
	int pNum;
	int cNum;
	int solverType;

#ifdef IMKL
	SVD_param svd;
#endif
};

#endif


